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FOREWORD 


This report describes the work performed by Business and Technological 
Systems, Inc. under Contract No. NAS 5-20045 with the NASA/Goddard Space Flight 
Center. The work deals with the development of a computer program to determine 
the ocean tidal perturbations in the elements of artificial earth satellites. 


1. INTRODUCTION 


Tidal forces have special significance in geophysics since their study can 
give insight into the properties of the earth's interior and its elastic re- 
sponses and ocean dynamics . Until recently the analysis of tidal perturbations 
of artificial satellites due to the jcftbined effect of earth and ocean tides 
has been applied to the determination of values for the characteristic parameters 
(Love numbers) which describe the elastic response of the earth with only the 
influence of solid earth tides taken into consideration. These results yield 
a value for the Love number k 2 somewhat smaller than anticipated and show an 
apparent dependence of k 2 on the orbital inclination. The perturbations due 
to ocean tides could be the primary cause of this disagreement. 

As the ocean tide, which manifests itself directly through gravitational 
attraction as well as indirectly through effects such as ocean loading is obser- 
vable in satellite data, a basic understanding of the mechanism Involved and 
an accurate ocean tidal theory is required for precise analysis of satellite 
perturbations based on the models of the principal ocean tidal components. 
Existing models of tidal component amplitudes and phase lags can be used to 
give a general estimate of oceanic effects on satellites. However, the compre- 
hensive treatment of the probelm requires the solution of the Laplace Dynamical 
Tidal Equation with realistic boundary conditions, friction and ocean loading 
to obtain more accurate tidal component models. 

In the present work we follow the semi-analytical treatment of Musen [5] 
for the oceanic perturbation equations and develop computer programs to calculate 
the influence of ocean tides on the satellite orbital elements in terms of the 
tidal amplitudes and lags of existing models for the M 2 and ocean tidal 
components. It is with pleasure that the author acknowledges the direction and 
support of Dr. Peter Musen of NASA/Goddard Space Flight Center in this effort. 


2. TIDAL COKSTITUENTb 



The sum of the gravitational potentials of the moon and sun at any point 
on the surface of the earth are periodic functions of time. In 1921 Doodson [2] 
presented a harmonic expansion based on the Brown lunar theory which Is the 
basis of modern tidal analysis. In particular, Doodson represented the tidal 
potential as a function of six basic astronomical variables in the form 

l A(ap . . . , a^) expi [a^t 4 a 2 s + a^h 4 a^p 4 a^N' 4 a & p g ] (1) 

a ^. • • • » d 2 

where 

r is mean lunar time 
s is the mean longitude of the moon 
h is the mean longitude of the sun 
p is the longitude of the lunar perigee 

N" is the negative of the longitude of the ascending lunar node 

p is the longitude of perihelion 
6 

and where the sum is taken over all positive and negative integer values of 
the a^. Each term in the expansion displaces the geoid a distance 

A(a } 

Sh - ~ 

g 

where g is the value of the local gravity, and these displacements are called 
tides or constituents. 

The tides fall into three major species denoted by a^ = 0 (long period), 

- 1 (diurnal) and = 2 (semidiurnal) , and are usually denoted by the 
Doodson number , defined as 

(a^, a 2 45, a 3 + 5, a^4 5, a 5 45, 8^45). 

Certain major constituents are labeled by letters, with a subscript denoting 
the integer a^. For example M 2> P^, K^, 0^ denote (2, 5, 5, 5, 5, 5), (1,6, 3, 5,5,5) , 
(1,6, 5, 5, 5, 5) and (1,4, 5, 5, 5, 5) respectively. 

The astronomical variables are very nearly linear functions of time and are 
expressed in mean solar days (T) from January 0, 1900 mean noon of Greenwich as 
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s * 270°. 43659 + 13° .1763967 T 
h « 279* . 69668 + 0?9856473 X 
p - 334*. 32956 + 0*. 11140408 T 
IT « 100“ .81G72 + 0 ,*05295391 
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281°. 22083 + 4! 70642026*10“ 5 T. 


Mean Lunar time (t) and mean solar time (t) are related to sidereal time 
(8) through the relations 

8 a t + h 
0 e x + 8 

so that 

x a t + h - s. 

and the variables A, l' t F, D, X of the lunar theory are related to the Doodson 
variables such that 

1 ® s - p 

r = h-p s 
f = s + rr 
D = s - h 

r " Ps 
or 

s = a' + d+ r 

h « A + r 

P - -nr + D + r 

= -r+ f - d - r 

p = r . 

*s 


For example, the M 2 tide (2, 5, 5, 5, 5, 5) has an argument 

2x + 0 • s + 0 • h + 0 • p + 0 • N ' + 0 • p = 2t + 2h ~ 2t - 2D 

s 

and the tide (1,6, 5, 5, 5,5) has 


t + s = 0 = t + h 


Typically the solutions for a tidal constituent a^=m i*_ represented in 
the form 

F m (x*<t>) cos (at + v + S(x,<J>)) (2) 

where the phase angle S(x,i)>) is the retardation of the high water at (x, <t>) 
relative to high water at Greenwich and t is the number of mean solar days from 

epoch t n . Tidal chart data is generally in the form of co-range lines (contours 

" f 7 1 

of equal amplitude) and co-tidal lines (isochrones of the retardation) 1 i with 

the time origin at lunar transit at Greenwich. 


1 


3. OCEAN TIDAL PERTURBATIONS 

As developed by Kusen^, the satellite perturbations due to ocean tides 
are expanded Into Fourier series with the orbital elements ft and 'S r of the 
satellite and 1, t', F, D and V of the lunar theory as arguments. The 
coefficients of the expansion are numerical and depend upon the satellite 
elements. For the sake of completeness we present a brief description of the 
Musen development. 

The disturbing function V can be written 

„ ™ 00 v n+1 . n+1 n 

V-3f f l <§) 


n=0 


u-Ht;) (?) 


y N P (sin 5) 
L n nm nm v 
ra^O 


(3) 


4 tr 


J ^(x - *, <T, t) P nm (sin <p') cos(m[a - 9 - x'])do' 


where G is the universal gravitational constant, M is the mass of the earth, 
r is the geocentric distance of the satellite, a is the semi— major axis of 
the satellite orbit, R is the radius of the mean level surface of the ocean, 

R£ is the height of the water tide, < 0 and k are the densities of the earth 
(mean) and sea water respectively, k' are Love numbers associated with the 
effects of loading of the tidal water mass, a and 5 are the right ascension 
and declination of the satellite, 8 is the sidereal time at Greenwich, 
are the associated Legendre polynomials and 


where 5 


ij 


j, _ 2 (n-m) 1 _ g 

nm (n+m) l m,0 


is the familiar Kronecker delta. 


(4) 


The tidal oscillations of the sea at a given point (x'» tp - *) on the mean 
sea surface is represented as a sum of periodic components of the form 


£(x% 4 , t) = l f = l F m (x“> <T) cos [m(8+x') +q m (x •',!()■') + 
♦ m m m 


(5) 


where ip is expressed as a linear combination of the lunar elements £, l , F, 
in 
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D and T, and ra denotes the distinct tidal component 0^, P^ f M 2 , etc.). 

This is to be compared with Eqn. (2) . The local sidereal time 0+x # is averaged 
out when only the long period effects are considered and the slowly changing 
part q + i|> m remains. 

The disturbing function then be written as a sum of components 

V - l V (6) 


where for a given constituent 


y B 3/2 — JL Y y 
V m ' R J* nm 
0 n'na 


R n+1 a n+1 

V = (— ) (l+k')A N Or) P (sin 6) cos (ma + ip ) 

mn a n nm nm r nm nm 


ti ip + e 

y nm - r m nm 


A nm cos e nm “ 47 F m (x% ^ )P nm C8in *' )c0S 


A m. 8ln 8 nm * 4? / j *,«<*'• ♦' )1 W 8ln ' n8ln • 


The long period effects are obtained from perturbation differential equations 
in which the short period terms are averaged over the instantaneous satellite 
orbit. Denoting 


Y__ * P _(sin <5) cos(ma+ i|» ) 

nm nm nm 


Z nm * tsin 1 sln a Jx 


rr - sin i cos ft ~ + cos i -^-]Y 


3v nm 


where X, p, v are the equatorial components of the geocentric unit vector in 


- 6 - 


1 


the direction of the satellite, we Itave 


2c 


d6e 


nm 


dt 


u/l-e 2 R 1 [ >a, P+1 ^ Y nm M 
B • (-) 1 OH 


e nm 2ir j r 
0 


3c r 


(13) 


2ir 


d$i 


nm 


n 


dt 


™ 2tr 


n+] 

,& " Z cos(£ +cr- 8)dK 
(“) nm 


(14) 


dsn 


n 


nm 


dt 


B 


/I^sini tUB 2 * 


2c 

f (£) lri ’ lz nm 8l n(f +^-»)dM 


(15) 


2 it 


d6ar_ 


nm 


dt 


» 1 
B 


/: — t nm 2 it 


n+1 


_ n-rx i , -» 

(|) / (n+1) (i e+cos f + ~e cos 2f)Y 


nm 


(16) 


3Y I d6fi 

+ (2 sin f+ ^ e sin 2f) > dM + 2 sin 2 i/2 1)31 


dt 


2c 


dSM 


nm 


dt 


2nB (n+1) * 
nm 


2ir 


n+1 d<5w 

(?) 


dt 


(17) 


+ 2/l-e^ sin 2 i/2 


d5fl 


nm 


dt 
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where 


B 


nra 


— ' c~) n (i+k:)A n n 
<q a n nm run 


(18) 


and where f is the true anomaly, e, nr, 0, a and i are mean elliptic elements 
of the satellite, n is the mean motion and H is the mean anomaly. 

The closed form integration of the perturbation equation*? (13) - (17) is 
facilitated by use of special purpose multi-dimensional Fourier series manipu- 
lation subroutines. For a specified tidal component m, Y nm and Z^ are 
expanded into Fourier series with arguments f, ter, 0 and iji to arbitrary 
values of n by recursion on functions of lower order and 1, y, v through 
the formulas 

Y n+l,m " n-m+1 ^? 2n+1 ? vY nm ~ Y n-l,m 

W ■ {< 2 " +1 >' ,Z n« + < 2 ” +1 >« 8 1 l n»- < 2 °> 

where 



X “ cos 2 i/ 2 cos (f+o) + sin 2 i/2 cos (f ■Fur- 2ft) 

y =» cos 2 i/2 sin(f+tzr) - sin 2 i/2 sin(f +«■- 2f!) (21) 

v ** sin i sin(fiur-n) . 

7ne coefficients in the expansions are then numerical and dependent on the 
satellite Inclination. 

The recursion is initiated by formally calculating 

J 

lor mF=l and 


Y^ = X cos - y sin tp 
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= sin i sin(fi + tfi) 


( 22 ) 


( 23 ) 


Y 22 “ 3(*‘ £ -U 2 ) c08 - 6Xw sin 


^22 * ^ 8 ^ u ^ sin O'i + '♦<)+ 6y sin 1 cos (& + <(/) 


for m-2 while noting that X^, 2 Jk are zero for j < k in Eqn. (iy). After 
all expansions are calculated is replaced by the proper value of 

nm 


The factors multiplying (^) in the integrands of Equations (13) - (17) 
are now finite Fourier series with arguments f, -ctr, ft and ^ . They are 
formally integrated term by term through the relation 


n+1 


1 ( a n+ l 

2 if J ( 7 } { sin } (Hff + Vr+V + ^dM 

(24) 

- d-« ! )' n+4 a^' 1 (a) • ('™> (n^+n^ + n^) 


v..»'re 



t 2 ?) 

I 

1=0 


, m . k+2i 
k+2 )e 


(k+21) 1 

(21) 1 J (2k+2£) ! 1 * 


(25) 


The resultant differential equations may now be Integrated in closed form to 
obtain the perturbations, since the arguments '■nr, ft and ij> of the averaged 
Fourier series are very nearly linear functions of time where 


and (Brouwer, 1959) 


n -n 0 + n(t-t 0 ) 

'tr «'nr 0 +nr(t- t Q ) 


. J. R cos i 

ft = - -in— - 

2 a 2 (1-e 2 ) 2 
J R 2 

. 3 J 2 R 2 

* - - fn-~ =-r 11+2 cos i-5 cos^ i] 

4 a Z (l-eV 


( 26 ) 
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4. PROGRAM FORMULATION 


A set of PL/I programs have been developed for the manipulation of the 
multi-dimensional Fourier series required for the perturbation equations. The 
storing, processing, combining, multiplying, integrating, averaging and evalu- 
ating of the series has been subrout ini zed and the main program performs the 
numerical calculations by a sequence of calls to these routines. The user has 
control of the program core size through variably dimensional arrays represent- 
ing the trigonometric series. The user specifies the minimum value of the 
series coefficients to be retained and NKAX, the maximum value for n, where 

n B m, mil, ... NMAX 

denotes the number of terms calculated in the expansion of the disturbing 
function, Eqr. (7). 

The scheme for storing the trigonometric series is as follows: Let y^ 

denote a trigonometric series associated with the calculation of perturbations 
for the n 1 term in the expansion of the disturbing function for constituent 
m such that 

SlAX 

Y n* * \ ® <V + "2 TO+ "3 a + n 4* ) - 

Then an integer array dimensioned 6) stores the trigonometric function 

and arguments while an array of dimension stores the coefficients. For 

a particular value of A, the first element of the integer array stores 0 
for a cosine term and 1 for a sine term, while elements 2 through 6 store 
n l* n 2’ n 3’ n 4 ant * n res P ectivei y- nu ®k er °f terms in the series, 

is stored separately as an integer variable. 

The program structure consists of a main program TIDES in which the user 
specif ie° the input for a particular satellite and for the tidal constituent 
desired, and the following subprograms: 

EXCUT: Calculates the ocean tidal perturbation on the satellite by 

recursively calculating Y nm and Z^, forming the perturbation 
differential equations, symbolically averaging the differential 
equations over the satellite orbit, integrating the equations in 
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II 


In closed form, and calculating for perturbations in elements 
over a given time interval. 


PAF: 


ANALIN: 


Computes the partial derivative of a trigonometric series with 
respect to f , ur or ft. 

n+1 

Calculates the average of a trigonometric series (“) y 

r nm 

over the instantaneous orbit of the satellite. 


C0MBAL: Computes the combinatorial coefficient (“) . 


SERMUL: Multiplies two trigonometric series together. 

SERSHM: Adds two trigonometric series together. 

HANSEN: Calculates the Hansen coefficient a^(e) by recursion. 


EVAL1: Prints a trigonometric series and the period of each term. 

LUNARG: Computes the arguments , F, D, r at time t. 

EVALDS: Evaluates a series at time t. 

COMP: Compresses a trigonometric series by deleting terms whose coefficients 

are less than a specified tolerance. 

ENTSER: Integrates a series in closed form with respect to time. 

PASSER: A F0RTRAN subroutine which generates plots of the perturbations 

in elements. 

The parameters A , e of Equation (10) have been obtained from cotidal 
• nm . j .nm rgi 

and corange data for M^ L J and 1 from charts and included in the program. 
These charts are presented in Figures I and II. The chart contour data was 


obtained on magnetic tape from an electronic digitizing machine and then inter- 
polated onto a 3**3° global grid. The function q B (x, $) was calculated at the 
grid points from the cotidal data and the quadrature of Equation (10) over the 
global ocean was performed using a recursive two-dimensional integration scheme 
and a world function map which is zero on land masses and unity on the oceans. 
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5. PROGRAM OPERATION 


The PL/I program, as currently developed, operates for tidal constituents 

M 2 and K^. The investigation of other tides would required the input of the 

parameters A and t for those constituents in the program EXCUT. User 
nm nm ° 

input to the program is through the main program TIDES , where the initial 
conditions 


SI ® i (inclination in radians) ; 

ECC * e (eccentricity’ 1 , 

PIEO « rsf (longitude of perigee in radians); 

0MEGO * 0 (longitude of ascending node in radians); 

L0NG ® L ** to' 1* M (mean longitude in radians); 

SN >- n (mean motion in radians/day) ; 

P = 7 (reciprocal semi-major axis in earth radii); 

a 

TJD = Julian date at epoch; 

are specified for a particular satellite. The tidal constituent is defined by 


M 


2 £orM 2 "\ _ 

1 for K x | 


while the constituent argument = n^S. + + n^D + n^T is given by 
the integer array NPSI(5) where 


NPi’I(i) => n ± ; 


and \p - 
m 


from 


~2l'-2D-2F for M_ and tj/ =0 for K-, . 

z m 1 

The recursion for the terms of the disturbing function is continued 
n = m, m+1, NMAX 


and terms in the perturbation series are deleted for values of the coefficients 
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less than T0LER. Both NMAX and T0LER are user specified. Additionally, the 
integer MAXDM defines the dimensionality of the arrays storing the trigonometric 
series and is a function of NMAX and T0LER. As NMAX or T0LER (or both) become 
larger, more terms will appear in the perturbation series so that MAXDM must 
be increased* 

The variables ASTART, AFIN and A1NC specify the interval in days 
about the epoch TJD for which perturbations will be calculated. This interval 
is (TJD + ASTART, TJD + AFIN) in increments of A1NC days. Note that these varia- 
bles are not necessarily Integers and may be positive or negative. 


6. NEW TECHNOLOGY 


The effort under this contract consisted of the development and programming 
of techniques to calculate the Influence of ocean tides on satellite orbital 
elements in terms of the tidal amplitudes and lags of existing constituent 
models. High order expansion of the oceanic tidal potential was obtained by 
recursion with multi-dimensional Fourier series and the perturbation differen- 
tial equations averaged over the satellite orbit and analytically integrated. 
The parameters needed for the and tides were obtained from global inte- 
gration of tidal data taken from existing charts. The and perturbations 
acting on the BE-C satellite were evaluated over a one hundred day interval 
and compared with the perturbations due to solid earth tides. 

Frequent reviews and a final survey for new technology were performed. It 
is believed that the mathematical and programming techniques and algorithms 
developed do not represent "reportable items," or patentable items, within the 
meaning of the New Technology Clause. Our reviews and final survey found no 
other items which could be considered reportable items under the New Technology 
Clause. 



FiGL-ac 10. Cotidat lines ( ) andcorangc lines (metres) ( ) for the Af, tide; map 1* A; computation grid, Jt = l*;/= <rj?t j = O.OC33; F^ptr.h^ ~ U 

(a = 0.5, A, = 1 km, n = 2, F u = frictional foicc/unit area). 






7. 


NUMERICAL COMPUTATIONS 


The and tidal perturbation equations were expanded to order n*20 
for the BE-C satellite with the following initial conitlona: 

e - .025037 

i - 41U929 

a ■ 1.177 (earth radii) 

n = 301! 712 

/nr « 272!139 

M - 353!l64 

n “ 83.8567 (radians per day) 

Tq = 2440812.5 (Julian days). 

These expansions produced no terms in the expressions for the perturbed elements 
for n greater than 14 whose coefficients were within two orders of magnitude 
of the principal terms of the perturbation. Consequently, for the results pre- 
sented in this report the expansions were truncated at n=15. The principal 
terms for the perturbed elements are listed in Tables II and III, while Figures 
III - VII display the M 2 and K x perturbations of the BE-C satellite orbit 
evaluated over a one hundred day interval compared to the perturbations due to 
solid earth tides (k 2 = .29) ^ . 

The M 2 and tides produce an appreciable affect on the orbital perturbations. 
In the inclination these perturbations can be approximately 10 per cent of those 
due to solid earth tides. Although the M 2 tide has a larger amplitude than 
over portions of the world ocean, has a greater influence on the orbit of 
BE-C. Computations for the BE-C orbit has shown that the indirect effect of 
the M„ and K. ocean tides with J 0 has a negligible influence, in contrast to 

*■ rgi *■ 

the solid earth tides. J 

It must be noted that these results are of a preliminary nature and are 
dependent on the perturbation parameters A^ and calculated from existing 
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chart data. These calculations require the co-range and co-tidal values at 
global grid points and thus demand extrapolation in some regions of sparse 
chart contour data. In addition, to more fully understand the total effects 
of ocean perturbations acting on satellites, other influential constituents 
such as 0^, P^, S 2 , etc., must be investigated. This suggests the need for a 
new accurate numerical integration of the Laplace Tidal Equations with realistic 
force models for the important tides on a uniform basis. The knowledge of the 
ocean perturbations due to the superposition of the major tidal constituents 
will then lead to improved values of the earth's elastic parameters (Love 
numbers) from satellite observations. 
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n 

A * R_ 
nm & 

M 2 

(centimeters) 

e (radians) 
nm 

A * K 

nm l 

K 1 

(centimeters) 

t (radians) 

nm 

1 




.7894 


1.069 

2 

25.629 


-1.89 

2.583 


0.008 

3 

11.556 


-0.042 

2.724 


2.357 

4 

55.309 


1.308 

2.509 


2.609 

5 

55.497 


-2.282 

1.580 


-1.099 

6 

13.856 


-3.05 

0.760 


-0.567 

7 

45.265 


-1.649 

2.164 


-1.015 

8 

22.731 


1.849 

0.145 


2.573 

9 

29.368 


-0.615 

0.847 


3.113 

10 

17.565 


1.148 

0.476 


2.563 


TABLE I. Perturbation parameters obtained from the 

M2 (Pe ,r eris and Accad) and (Zahel) tides. 
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Perturbation 


Terms 


Period (days) 


Se 

-1.836*10 8 sin (iff + U + 4*22^ 

12.13 


+2.439*10" 8 sinenf- 3«- v 52 ) 

8.99 


+1.093*10“® sinftff+ii + tiO 

12.13 

6i 

1.2734*10" 7 cos(2s) + vf/ 22 ) 

10.33 


- .9254*10 7 cos(2Q+ ij/^ 2 ) 

10.33 

6 n 

-1.455*10" 7 8in(20 + «i> 22 ) 

10.33 


+ .1489*10" 7 sln(2R + tp c .) 

10.33 

finr 

-7, 339*10" 7 cos('nr+n + ^ 32 ) 

12.13 


-1.992*10" 7 cos (iff- 3n-* 32 ) 

8.99 


+9.760*10 7 cos (or- 3fi - ii^ 2 ) 

8.99 


+4.382*10 7 cos (ttf+ + ip 72 ) 

12.13 


-2. 761*10" 7 cos(lf-3fl-i|» 72 ) 

8.99 


-1.252*10“ 7 cos(-tff+fi + ^ 11 2 ) 

12.13 


+1.027*10~ 7 cos (nr- 30 - if> u 2 ) 

8.99 

6L * 5nJ + 6M 

-2. 875*10~ 7 sin(2fl + t(i 22 ) 

10.33 


+3.040*10" 7 sin(2fi + i|j 42 ) 

10.33 


Table II. Principal terms from the expansion for perturbed elements due to 
the M 2 tide for the BE-C satellite. 
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Perturbation 

Terms 

Period (days) 

fie 

3.497*10" 8 cosCw- 2n-« 31 ) 

38.27 


-1. 395*10“ 8 cos(u)-2n-t 51 ) 

38.27 


-2.144*10 8 cos(u + ifiy^) 

391.38 


-2.803*10 8 .cos (w 4* 

391.38 


+1.113*10“ 8 cos(n + iP 1;L ± ) 

391.38 

fii 

-2. 408*10" 7 sin(n + ij>2^) 

84.84 


3.843*10~ 8 sin(n + »i/ 41 ) 

84.84 

fin 

-6.439*10 8 eosCQ + ^j) 

84.84 


-2. 664*10“ 7 cos(n + ^ 41 ) 

84.84 

5m* 

-1.398*10 8 sin (w - 2 n - 

38.27 


—6 

-6.060*10 ain(w + tjj^) 

391.38 


+8. 594*10" 7 sin(u + ^ 71 ) 

391.38 


+1.126*10“ 6 sin(fi + ^ gl ) 

391.38 

5L » 6ttr'+ fiM 

-9.674*10 7 cosCn + tl^) 

84.84 


Table XXX. Principal terms from the expansion for perturbed ilements due 
to the tide for the BE-C satellite. 
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